/*
 * Main program
 * Copyright (C) 2010 Vsevolod Scherbinin
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <netcdf.h>
#include <assert.h>
#ifndef ADT_COMPLEX
#include "adt_complex.h"
#define ADT_COMPLEX
#endif

//#define MAX(a, b) ((a) > (b) ? (a) : (b))
//#define MIN(a, b) ((a) < (b) ? (a) : (b))

#define ADTFUNTYPE
#define ADTFUNNAME(fn) adt_complex_##fn
#define ADTFUNNAMEB(fn) adt_bicomplex_##fn
//#ifndef ADT_COMPLEX_INTERNAL
#include "adt_complex_internal.c"
//#define ADT_COMPLEX_INTERNAL
//#endif
#undef ADTFUNTYPE
#undef ADTFUNNAME
#undef ADTFUNNAMEB

#define ADTFUNTYPE __device__
#define ADTFUNNAME(fn) adt_cuda_complex_##fn
#define ADTFUNNAMEB(fn) adt_cuda_bicomplex_##fn
#include "adt_complex_internal.c"
#undef ADTFUNTYPE
#undef ADTFUNNAME
#undef ADTFUNNAMEB

#define CUBLOCKSIZE 64

#define ADTFUNNAME(fn) adt_cuda_complex_##fn

//This function is intended for error message generation. 
//Copying from original documentation
void checkCUDAError(const char *msg)
{
  cudaError_t err = cudaGetLastError();
  if( cudaSuccess != err) 
  {
    fprintf(stderr, "Cuda error: %s: %s.\n", msg, 
	    cudaGetErrorString( err) );
	    exit(EXIT_FAILURE);
  }                         
}

//Integrand for Hy
__device__ adt_complex_float
intgHy(float x, adt_complex_float alph, adt_complex_float w, float sinthet, float xi, float zeta)
{
  //alph - variable on contour; w - transversal wavenumber
  //den - integrand demominator; exx - integrand numerator
  adt_complex_float tmp, den, exx, cUnit, aux;
  ADT_SET_COMPLEX(&cUnit,1.,0.);
  //   den = (alph - sinthet)*csqrt(1. - alph);
  den = ADTFUNNAME(sqrt)(ADTFUNNAME(sub)(cUnit,alph));
  den = ADTFUNNAME(mul)(ADTFUNNAME(sub_real)(alph,sinthet),den);
  //   exx = cexp(6.28318530718*I*(fabs(zeta)*w - xi*alph));
  aux = ADTFUNNAME(mul_real)(alph,xi); //xi*alph
  tmp = ADTFUNNAME(mul_real)(w, fabs(zeta)); //fabs(zeta)*w
  tmp = ADTFUNNAME(sub)(tmp,aux); //fabs(zeta)*w - xi*alph)
  exx = ADTFUNNAME(exp)(ADTFUNNAME(mul_imag)(tmp,6.28318530718));
  return ADTFUNNAME(div)(exx,den);
}

// Kernel with integration
__global__ void
HalfPlane( adt_complex_float* cuHy, adt_complex_float* cuSf, float* cuxdt, float* cuzdt, float* d_xG, float* d_wG, float cukap, float cucos, float cusin, float Mult, int Xin, int NN, int Fract)
{
  int i = blockIdx.x*CUBLOCKSIZE+threadIdx.x; //Thread number defining
 int nodeindx, j;
 __shared__ float xi, zeta, Sx, Sz, Snrm;
//float d_xG[2000];
// float d_wG[2000];
  adt_complex_float Iw, tmp, auxH, auxX, auxZ, crHy, crSf, Ex, Ez;
  adt_complex_float cUnit, intHy, alph, w;
  
  ADT_SET_COMPLEX(&cUnit,1.,0.);
  
  if (i < NN) {
     //Coordinates X and Z picking up.
     xi = cuxdt[i%Xin]; zeta = cuzdt[i/Xin];
    
    //Incident wave computing
    tmp = ADTFUNNAME(rect)(6.28318530718*(xi*cusin + zeta*cucos),0.);
    Iw = ADTFUNNAME(exp)(ADTFUNNAME(mul_imag)(tmp,-1.));
    
    //  cHy = Iw - I*mult*sign(zeta)*complGauss16(-hlim, hlim, fract, intgHy); 
    ADT_SET_COMPLEX(&auxH,0.,0.);
    ADT_SET_COMPLEX(&auxX,0.,0.);
    ADT_SET_COMPLEX(&auxZ,0.,0.);
    nodeindx = 16*Fract;
    for (j=0; j<nodeindx; j++) {
      //   alph = x + kappa*I; w = csqrt(1. - alph*alph);
      //   w = creal(w) + cabs(cimag(w))*I;
      ADT_SET_COMPLEX(&alph,d_xG[j],cukap);
      w = ADTFUNNAME(sqrt)(ADTFUNNAME(sub)(cUnit,ADTFUNNAME(mul)(alph, alph)));
      ADT_SET_IMAG(&w,fabsf(ADT_IMAG(w)));
      intHy = intgHy(d_xG[j], alph, w, cusin, xi, zeta);
      auxH = ADTFUNNAME(add)(auxH,ADTFUNNAME(mul_real)(intHy,d_wG[j]));
      auxX = ADTFUNNAME(add)(auxX,ADTFUNNAME(mul_real)(ADTFUNNAME(mul)(w,intHy),d_wG[j]));
      auxZ = ADTFUNNAME(add)(auxZ,ADTFUNNAME(mul_real)(ADTFUNNAME(mul)(alph,intHy),d_wG[j]));
    }
    
    crHy = ADTFUNNAME(sub)(Iw,ADTFUNNAME(mul_imag)(auxH,Mult*ADT_SGN(zeta))); // Iw - tmp*I*mult*sign(zeta)
    ADT_SET_COMPLEX(&crSf,0.,0.);
    
    Ex = ADTFUNNAME(add)(ADTFUNNAME(mul_real)(Iw, cucos), ADTFUNNAME(mul_imag)(auxX,Mult));
    Sz = ADT_REAL(ADTFUNNAME(mul)(Ex,ADTFUNNAME(conjugate)(crHy)));
    Ez = ADTFUNNAME(add)(ADTFUNNAME(mul_real)(Iw,-cusin), ADTFUNNAME(mul_imag)(auxZ,Mult*ADT_SGN(zeta)));
    Sx = ADT_REAL(ADTFUNNAME(mul)(Ez,ADTFUNNAME(conjugate)(crHy)));
    
    Snrm = sqrtf(Sx*Sx+Sz*Sz); Sx = -Sx/Snrm; Sz = -Sz/Snrm;
    ADT_SET_COMPLEX(&crSf,Sx,Sz);

    cuHy[i] = crHy; cuSf[i] = crSf;
 }
  
  __syncthreads();
}  

#undef ADTFUNNAME

// Host code
void gaussnodes(float* nodes, float* weights, int nn, float ll, float hh) {
  float wg[8] = {0.027152459411754094852, 0.062253523938647892863, 0.095158511682492784810, 0.124628971255533872052, 0.149595988816576732081, 0.169156519395002538189, 0.182603415044923588867, 0.189450610455068496285};
  float xg[8] = {0.989400934991649923596, 0.944575023073232576078, 0.865631202387831743880, 0.755404408355003033895, 0.617876244402643748447, 0.458016777657227386342, 0.281603550779258913230, 0.095012509837637440185};
  float a, b, step, mvr, nrm;
  int i, k, imv;
  step = (hh-ll)/(float)nn;
  a = ll;
  for (i=0; i<nn; i++){
    b = ll + step*(float)(i+1);
    mvr = 0.5*(b + a);
    nrm = 0.5*(b - a);
    imv = i*16;
    for (k=0; k<=15; k=k+2) {
      nodes[imv+k] = mvr + nrm*xg[k];
      nodes[imv+k+1] = mvr - nrm*xg[k];
    }
    for (k=0; k<=7; k=k++) {
      weights[imv+k] = wg[k]*nrm;
      weights[imv+15-k] = weights[imv+k];
    }
    a = b;
  }
}

#define ADTFUNNAME(fn) adt_complex_##fn

static void
complex_abs_arg (float *abs, float *arg,
		 const adt_complex_float *ary, size_t elts)
{
  size_t rest;
  float *mp, *ap;
  const adt_complex_float *sp;
  for (rest = elts, mp = abs, ap = arg, sp = ary;
       rest > 0;
       rest--,
	 *(mp++) = ADTFUNNAME(abs) (*sp),
	 *(ap++) = ADTFUNNAME(arg) (*sp),
	 sp++)
       ;
}

static void
complex_re_im (float *abs, float *arg,
		 const adt_complex_float *ary, size_t elts)
{
  size_t rest;
  float *mp, *ap;
  const adt_complex_float *sp;
  for (rest = elts, mp = abs, ap = arg, sp = ary;
       rest > 0;
       rest--,
	 *(mp++) = ADT_REAL(*sp),
	 *(ap++) = ADT_IMAG(*sp),
	 sp++)
       ;
}

int
main (int argc, char * argv[])
{

  FILE *logfile;
  FILE *cfgfile;

  // Input and output variables
  int fract, xin, zin, xzin; //Xin, Zin - Numbers of points for field computing
  float xmax, zmax, theta; //theta - incident angle (degrees) counted out z axis; Xmax, Zmax - Maximum values of X and Z in wavelength
  //number of fractions for integration
  float kappa, mult, hlim;
  // Internal variables
//   int narg; //Number of command line arguments
  int nb_integral, nb_xdt, nb_zdt, nb_fields; //Number of bytes for integration nodes and weights, for X points, for Z points and for fields.
  int xcnt, zcnt; //Counters for coordinates
  
  char basename[40], configfile[40], ncname[40];
  char str[50], prefix[6];

  adt_complex_float * Hy;
  adt_complex_float * Sf;
  float * xdts;
  float * zdts;
  float * h_xG;
  float * h_wG;
  float xis, zis, sinthet, costhet;
  float gpuTime = 0.0f;
  
  time_t timestart;
  clock_t timess, timeff;
  double cpuTime;
  

//===================================================================
//	COMMAND LINE ARGUMENTS PROCESSING
//===================================================================
  strcpy(configfile, "halfplane_compar.ini");
  if(argc>1){
    strcpy(configfile, argv[1]);
  }
//===================================================================
//	FILES OPENING AND VARIABLES INITIALIZING
//===================================================================
  xin = 641; zin = 640;
  xmax = 4.; zmax = 4.;
  
  logfile = fopen("current.log", "r");
  cfgfile = fopen(configfile, "r");


start:
	  while(fgets(str, 100, cfgfile)) {
    if ((str[0]=='s')&&(str[1]=='t')&&(str[2]=='o')&&(str[3]=='p')) {
      printf ("%s\n", "Stop found!");
      goto finish;
    }
    if ((str[0]==';')||(str[0]=='#')) {
//       printf ("%s\n", "Comment found!!!!");
      goto start;
    }
    time(&timestart);
    timess = clock();

    strcpy(basename, str);
    if(basename[strlen(basename)-1]=='\n') basename[strlen(basename)-1] = '\0';
    
    printf ("\nTask %s starts at %s", basename, ctime(&timestart));
    
    {
    int devID;	
    cudaDeviceProp props;
    cudaError_t cuderr;

    //Get GPU information
    cudaGetDevice(&devID);
    cudaGetDeviceProperties(&props, devID);
    printf("\tDevice %d: \"%s\" with Compute %d.%d capability; ", 
	  devID, props.name, props.major, props.minor);
    
    }
    
    //Data reading
    fscanf(cfgfile, "%f\n", &theta);
    fscanf(cfgfile, "%f\n", &xmax);
    fscanf(cfgfile, "%f\n", &zmax);
    fscanf(cfgfile, "%d\n", &xin);
    fscanf(cfgfile, "%d\n", &zin);

    //=================================================================
    //	Data preparing
    //=================================================================
    //Xmax and Zmax must bee positive
    xmax = fabs(xmax); zmax = fabs(zmax);
    hlim = 4.*MAX(xmax,zmax);
   fract = ((int)(2.*hlim))*10;
	//Number of fraction for numerical integrating
  theta = theta*0.0174532925199;	//Incident angle to radian
    sinthet = sinf(theta);
    costhet = cosf(theta);
    //Multiplyer 1D-1 cannot be minimized, cause pole set on integration contour
    kappa = (sinthet - 1.)*1E-1;
    //Zin must be positive, even and not less than 2
    zin = 2*(int)(abs(zin)/2);
    if (zin<=2) {
      zin = 2;
    }
    //Xin must positive, odd and not less than 3
    xin = 2*(int)(abs(xin)/2)+1;
    if (xin<=3) {
      xin = 3;
    }
    //File names preparing. Prefix needs to determine area for plotting.
    //Files for ampl[itude], phas[e] and [power ]flow.
    sprintf(prefix, "m%d_", (int)xmax);
    
    //Array of field and dots making
    xzin = xin*zin;
    nb_fields = xzin*sizeof(adt_complex_float);
    nb_integral = 16*fract*sizeof(float);
    nb_xdt = xin*sizeof(float);
    nb_zdt = zin*sizeof(float);
    
    Hy = (adt_complex_float *) malloc(nb_fields);
    Sf = (adt_complex_float *) malloc(nb_fields);

    xdts = (float *) malloc(nb_xdt);
    zdts = (float *) malloc(nb_zdt);
//cudaMallocHost( (void**)&Hy, nb_fields);
//cudaMallocHost( (void**)&Sf, nb_fields);
    //==============================================================================
    //	Computing
    //==============================================================================
    //Data preparing: computing normalizing multiplier and coordinate grid
    mult = costhet/(6.28318530718*sqrt(1.+sinthet));

    xis = xmax/(float)((xin-1)/2);
    for (xcnt = 0; xcnt <=xin-1; xcnt++){
      xdts[xcnt] = -xmax + xcnt*xis;
    }
    zis = zmax/(float)(zin/2);
    for (zcnt = 0; zcnt <= (zin/2-1); zcnt++) {
      zdts[zcnt] = -zmax + zcnt*zis;
    }
    for (zcnt = zin/2; zcnt<=zin-1; zcnt++) {
      zdts[zcnt] = -zmax + (zcnt+1)*zis;
    }

    h_wG = (float *) malloc(nb_integral);
    h_xG = (float *) malloc(nb_integral);
    gaussnodes(h_xG, h_wG, fract, -hlim, hlim);
    
    //==============================================================================
    //Field computing starts
    //device memory allocation
    {
    adt_complex_float * cuHy = NULL;
    adt_complex_float * cuSf = NULL;
    float * cuxdt = NULL;
    float * cuzdt = NULL;    
    float * d_xG = NULL;
    float * d_wG = NULL;
    
    // allocating GPU memory
   cudaMalloc ( (void**)&d_xG, nb_integral );
   cudaMalloc ( (void**)&d_wG, nb_integral );
    cudaMalloc ( (void**)&cuxdt, nb_xdt );
    cudaMalloc ( (void**)&cuzdt, nb_zdt );
    cudaMalloc ( (void**)&cuHy, nb_fields );
    cudaMalloc ( (void**)&cuSf, nb_fields );

    checkCUDAError("allocation");
    
    // create cuda event handles
    cudaEvent_t start, stop;
    cudaEventCreate ( &start );
    cudaEventCreate ( &stop );

    // asynchronously issue work to the GPU (all to stream 0)
    cudaEventRecord ( start, 0 );
    
    cudaMemcpy      ( cuxdt, xdts, nb_xdt, cudaMemcpyHostToDevice );
    cudaMemcpy      ( cuzdt, zdts, nb_zdt, cudaMemcpyHostToDevice );
   cudaMemcpy      ( d_xG, h_xG, nb_integral, cudaMemcpyHostToDevice );
  cudaMemcpy      ( d_wG, h_wG, nb_integral, cudaMemcpyHostToDevice );
    
    checkCUDAError("loading data");

    int threadsPerBlock = CUBLOCKSIZE;
    int blocksPerGrid = (xin*zin + threadsPerBlock - 1) / threadsPerBlock;
    printf("blocks %d or %d.\n", (xin*zin + threadsPerBlock - 1)/threadsPerBlock, 1 + ((xin*zin - 1)/threadsPerBlock));
    
    HalfPlane<<<blocksPerGrid, threadsPerBlock>>> ( cuHy, cuSf, cuxdt, cuzdt, d_xG, d_wG, kappa, costhet, sinthet, mult, xin, xin*zin, fract );

    
    cudaThreadSynchronize();
    checkCUDAError("syncing after calc");
    

    cudaMemcpy      ( Hy, cuHy, nb_fields, cudaMemcpyDeviceToHost );
    checkCUDAError("copy cuHy");
    cudaMemcpy      ( Sf, cuSf, nb_fields, cudaMemcpyDeviceToHost );
    checkCUDAError("copy cuSf");
    
    cudaEventRecord ( stop, 0 );

    cudaEventSynchronize ( stop );
    cudaEventElapsedTime ( &gpuTime, start, stop );

    cudaEventDestroy ( start  );
    cudaEventDestroy ( stop  );
    cudaFree         ( d_xG  );
   cudaFree         ( d_wG  );
   cudaFree         ( cuxdt  );
    cudaFree         ( cuzdt  );
    cudaFree         ( cuHy  );
    cudaFree         ( cuSf  );
    
    checkCUDAError("GPU memory clearing");
    
//     printf("Computing ends...\n");
    
    gpuTime*=0.001;
    }
    
    //
    //	The next block devoted to output the results to the NetCDF file
    //
    {
   int status;
    int ncout, xdim, zdim, absHy_id, phaHy_id, outSx_id, outSz_id, xpnts_id, zpnts_id;
    int ang_id;
    size_t sizz;
    int  nc_dimids[2];             
    
    snprintf(ncname, sizeof(ncname), "%s%s", basename,".nc");
    status = nc_create(ncname, NC_SHARE, &ncout);
    if (status != NC_NOERR) {
      printf("File %s exists.\n",ncname);
      status = nc_open(ncname, NC_WRITE, &ncout);
    }
    sizz = xin;
    status = nc_def_dim(ncout, "x", sizz, &xdim);

    printf ("\txin=%d; zin=%d\n", xin, zin);
    sizz = zin;
    status = nc_def_dim(ncout, "z", sizz, &zdim);
    nc_dimids[0] = zdim;
    nc_dimids[1] = xdim;

    status = nc_def_var(ncout, "xpoints", NC_FLOAT, 1, &xdim, &xpnts_id);
    status = nc_def_var(ncout, "zpoints", NC_FLOAT, 1, &zdim, &zpnts_id);
    status = nc_def_var(ncout, "absHy", NC_FLOAT, 2, nc_dimids, &absHy_id);
    status = nc_def_var(ncout, "phaHy", NC_FLOAT, 2, nc_dimids, &phaHy_id);
    status = nc_def_var(ncout, "flwSx", NC_FLOAT, 2, nc_dimids, &outSx_id);
    status = nc_def_var(ncout, "flwSz", NC_FLOAT, 2, nc_dimids, &outSz_id);

    status = nc_put_att_float(ncout, NC_GLOBAL, "inc_angle", NC_FLOAT, 1, &theta);
    status = nc_put_att_float(ncout, xpnts_id, "maximum", NC_FLOAT, 1, &xmax);
    status = nc_put_att_float(ncout, zpnts_id, "maximum", NC_FLOAT, 1, &zmax);
    
    status = nc_enddef(ncout);
    
    status = nc_put_var_float(ncout, xpnts_id, xdts);
    assert(status==0);
    status = nc_put_var_float(ncout, zpnts_id, zdts);
    assert(status==0);

    {
      float absHy[xzin], phaHy[xzin], oSx[xzin], oSz[xzin];
      complex_abs_arg (absHy, phaHy, Hy, xzin);
      status = nc_put_var_float(ncout, absHy_id, absHy);
      assert(status==0);
      status = nc_put_var_float(ncout, phaHy_id, phaHy);
      assert(status==0);
      complex_re_im (oSx, oSz, Sf, xzin);
      status = nc_put_var_float(ncout, outSx_id, oSx);
      assert(status==0);
      status = nc_put_var_float(ncout, outSz_id, oSz);
      assert(status==0);
    }
    
    status = nc_close(ncout);
    }
    //
    //	Memory clearing
    //
//    free(Hy);
    Hy = NULL;
 //   free(Sf);
   Sf = NULL;
    
   free(xdts);
    xdts = NULL;
   free(zdts);
    zdts = NULL;  
    timeff = clock();
    cpuTime = ((double) (timeff - timess)) / CLOCKS_PER_SEC;
    printf ("Task %s has been computed during %.3f on CPU, %.3f on GPU and %.3f overall.\n", basename, cpuTime, gpuTime, cpuTime+gpuTime);
 //   fprintf (logfile, "Task %s has been computed during %.3f on CPU, %.3f on GPU and %.3f overall.\n", basename, cpuTime, gpuTime, cpuTime+gpuTime);
	
    goto start;
    
}
finish :  printf("\nAll jobs completing!..\n\n");
	  fclose(cfgfile);
	  fclose(logfile);
	  return 0;
}

#undef ADTPREFIX
